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Abstract. Noneqiuilibrium dynamics of rotating droplets are studied by molecular dynamics 
simulations. Small deviations from the theoretical prediction are observed when the size of a 
droplet is small, and the deviations become smaller as the size of the droplet increases. The 
characteristic timescale of the deformation is observed, and we find (i) the deformation timescale 
is almost independent of the rotating velocity with for small frequency and (ii) the deformation 
timescale becomes shorter as temperature increases. A simple model is proposed to explain the 
deformation dynamics of droplets. 



1. Introduction 

Droplets are small drops of liquid which are commonly observed in our life. Many things can 
be considered as ensembles of liquid droplets, such as rain, clouds, and splays, etc. Droplets are 
not only familiar to us, but their behaviors are also important for engineering, such as spray 
combustion [1], developments of Inkjet printers [2] and electronic sputtering [3], etc. Behaviors 
of droplets are mainly governed by two kinds of forces, namely, the surface tension and the 
inertial force. The surface tension works as restoring force trying to keep shape of a droplet, 
while the inertial force usually tries to deform and destroy it. A simple example of such balance 
between two forces can be seen in rotating droplets where the centrifugal force plays the role of 
the inertial force. One of the fundamental studies of such rotating droplets was investigated by 
Plateau [4]. A water-droplet was put in alcohol with same density in order to mimic a gravity- 
less system, and then the container was rotated. The droplet was gradually compressed to flat 
shape as the rotating velocity increased and became unstable when the angular velocity exceeded 
some critical value. While the experimental technique is quite simple, it is difficult to investigate 
this phenomena analytically because of the influence from the ambient fluid. Recently, levitated 
droplets have attracted much researchers' interests. A droplet is levitated by some external 
force, such as electromagnetic force [5], and then a containerless system is achieved which is 
convenient for comparison with theoretical predictions. The properties of rotating droplets in 



Figure 1. A typical snapshot of a static droplet. The white lines denotes the simulation box 
with the linear size L = 40. The total number of particles A'^ = 16000 and the temperature 
T = 0.6. The red and gray particles denote a liquid droplet and vapor, respectively. 



steady-state are investigated theoretically, and numerical works were followed [6, 7]. Also, the 
gravityless experiments were performed in spacelab by Wang et al. [8]. They reported that 
while the axisymmetric shapes are well described by the theoretical prediction, the bifurcation 
point from axisymmetric to nonaxisymmetric shapes locates at a lower rotation velocity than 
the theoretical prediction. Despite of such past studies, dynamics of the deformation, especially 
nonlinear motions which cannot be described as simple oscillations, have not been clarified yet. 
Dynamical aspects of the deformation, such as characteristic timescale of the deformation or 
how the droplets separates into fragmentations, are very important both for theoretical interests 
and for applications. Additionally, the validity of the continuum treatments should be tested 
for small droplets, since a droplet consists of many particles, and consequently, the surface of it 
has finite thickness, while thickness is usually ignored in the continuum theory. In the present 
study, we investigate the dynamics of the rotating droplets using molecular dynamics (MD) 
simulations. Using MD, dynamics of droplets can be studied naturally taking account of its 
thickness. We first make a brief review of the continuum treatments for rotating droplets in 
the steady state, and compare it to our results after describing details of numerical methods. 
Deformation dynamics are observed for different conditions, and a simple model is proposed in 
order to explain the results. 

2. Form of Rotating Droplet 

The form of rotating droplets are solved exactly by Chandrasekhar [9] when its shape is uniaxial. 
The steady forms of rotating droplets are determined by the balance between the surface tension 
and the centrifugal force. Here, we make the brief review of the form of the rotating droplet 
following the manner by Chandrasekhar. 

Consider a droplet rotating around z-axis in a system without the gravity. In the steady 
state, the gradient of the internal pressure of the droplet should be balanced by the centrifugal 



force as 



2 dp 

pa; r = — , (1) 

where p is the density of the droplet, u is angular velocity, p is the internal pressure, and r is 
the distance from the rotating axis, respectively. Equation (1) directly leads to 

p{r)=PQ + ^puP'r'^, (2) 

with the pressure pq at the center of the droplet (r = 0). At the surface of the droplet, the 
pressure should satisfy the following relation 

p = 7divn, (3) 

with the surface tension 7 and the normal vector of the droplet surface n. From Eqs. (2) and 

(3), we have 

Pq + -^piiP'r^ = 7divn. (4) 

While shape of a rotating droplet varies from the sphere, the shape will be uniaxial when the 
rotating speed is small. Therefore, the form of the droplet can be expressed as z = /(r), where 
z is the height from the equatorial plane and r is distance from 2;-axis, namely, = a;^ + y^, 
respectively. Then the unit normal vector of the droplet surface n = {nx,ny,nz) is 
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Prom the above, we easily have, 
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Substituting Eq. (11) to Eq. (4), we obtain 
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which is integrable to be 
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Note that, the constant of the integration should be zero because of the rotational S3mimetry 
around z-axis. Let Re be the equatorial radius of the rotating droplet. Since (f) — >■ —oo when 
r ^ Re, we have 

27 87 ■ ^ ' 

We introduce a dimensionless value S as 



87 

where ojq is the characteristic frequency which is defined as 



Hereafter, we measure the length in the unit of the equatorial radius as f = r/Re- Equation (13) 
is then reduced to be 

^==f = -Kl-S + Sr-2), (17) 

or equivalently, 

= -^=^, (18) 

vl -r 

where 

ff(f) =f(l-S + Sf2). (19) 
Finally, we obtain the form of the rotating droplet by the following integration. 



f{f)=r<i)df. (20) 



While Eq. (20) can be expressed by the Jacobi elliptic functions, the expression is inconvenient 
for comparison with numerical results. We, therefore, derive the expression for the ratio of the 
equatorial radius of the rotating droplet to the radius of the droplet at rest. The volume of the 
droplet V is expressed in terms of /(f) by 

V = A'kRI ff{f)df. (21) 

The integration by parts leads to 
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= 27r['^L=df, (23) 
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= 27r/i, (24) 

where h is the function of the reduced rotational frequency uj/ujq. Consider the sphere with the 
radius Rq which volume is the same as the rotating droplet, that is, V = AttRq/S. Provided 
that the volume are not changed by rotation, Rq denotes the radius of the droplet at rest. The 
radius ratio Re/Ro is then expressed as, 

i?e/i?o=y . (25) 
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Figure 2. The density of the droplet for = 16000. The origin is set at the center of the 
droplet. The solid line denotes the hyperbolic tangent function of the form tanh((r — Rq)/X) 
with the radius of the droplet Rq = 15.7(3) and the thickness of the surface A = 0.88(5). The 
ratio of the thickness to the radius of the droplet is about 5.6%. 



Equation (25) means that deformation behaviors can be scaled for different rotation speed, 
volume of droplets, surface tension, and so forth. Equation (25) is shown as the solid line in 



A couple of things are worth to be noted. First, the equilibrium form of the rotating droplet 
can be given by the variational principle with appropriately chosen effective potential function, 
as described in Ref. [6]. An Euler-Lagrange equation derived by the variational principle leads 
to the Young-Laplace equation (3). Therefore, the above arguments are completely equivalent 
to those from the variational principle. Second, one has to be careful with the definition of the 
characteristic frequency defined in Eq. (15), since some researchers take the radius of the droplet 
at rest Rq as the scaling length, while we use the equatorial radius of the rotating droplet i?e- 
Accordingly, the definition of S and ujq can be different for researchers. 

3. Method 

In order to study the dynamics of rotating droplets, we perform MD simulations. We use the 
truncated Lennard-Jones potential of the form 



with the well depth e, the atomic diameter a, and the cut-off length rc [10]. The coefficients 
Co and C2 are determined so that V{rc) = V'{rc) = 0, i.e., the values of potential and the force 
become continuously zero at the truncation point. We choose the cutoff-length as rc = 3.0. In 
the following, we measure the physical quantities in the unit of the radius a, the well depth e, 
the Boltzmann constant k^, and the particle mass m. Time step is chosen to be 0.005. The 
simulation box is the cube with the linear size L, and the periodic boundary condition is taken for 
all directions. Total number of particles N = 4000, 8000 and 16000 are studied. Temperature is 
controlled by Langevin thermostat [11]. Simulations are partially performed by MDACP which 
is freely available online [12] 
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Table 1. Physical quantities of droplets at rest for N = 4000,8000, and 16000. Rq is the 

position of Gibbs surface, p is the density of the liquid, A is the thickness of the surface, and 7 
is the surface tension, respectively. While the density and the thickness are almost independent 
of the size of the droplets, the value of the surface tension increases as the size of the droplet 
increases. 

Particles are initially distributed spherically at the facc-centered-cubic lattice, and time 
evolution is performed under the fixed temperature T = 0.6. The density of the total system 
is p = N/L^ = 0.25, which is the liquid- vapor coexistent phase at the temperature. The 
system, therefore, reaches its equilibrium state after sufficiently long time, and the system 
contains a single droplet in the saturated vapor. In the present study, 10® steps are spent 
for the thermalization. An typical snapshot of a droplet is shown in Fig. 1. In order to identify 
the droplet, we define that two particles within the distance 1.3 belong to the same cluster, 
and the largest cluster is defined as the droplet which is shown in red particles in the following 
figures. The density profile is shown in Fig. 2. Since the interface between the liquid and the 
gas has finite thickness, we have to define the equatorial radius to compare the results with 
the theory. Here, we use the position of Gibbs surface to define the equatorial radius. Let r 
be distance from the center of a droplet. The local density p{r) is well approximated by the 
hyperbolic tangent function as, 

p(r) = (pL - Pg) tanh ((r - Rg)/X) + PG, (27) 

with the liquid density pL, the gas density pc, the position of Gibbs surface Rq, and the thickness 
of the surface A. From the fitting using Eq. (27), we determine the values of pL, PG; and Re- 
Then we identify the position of Gibbs surface Rq as the equatorial radius i?c. We measure 
the surface tension of the droplets in this equilibrium state following the method proposed by 
Ikeshoji et al. [13]. When the droplet is at rest, the equatorical radius Re equals Rq. The 
measured physical quantities of the droplet at rest arc summarized in Table. 1. 

After the equilibrium droplet is obtained, the thermostat is turned off and angular velocity 
u) is given to the droplet. As the droplet deforms, the moment of inertia also changes, and 
consequently, the angular velocity varies. It makes difficult to investigate the relation between 
the angular velocity and the final form of the droplet precisely. Therefore, we keep the angular 
velocity fixed throughout the time evolution by the method similar to the velocity scaling 
scheme [14]. We observe the angular velocity of the droplet, and rescale the velocities of the 
particles in the droplet so that the total angular velocity of the droplet is kept to be the initial 
value. 

4. Results 

4-1- Static Form of Rotating Droplets 

Two typical time evolutions of droplets are shown in Figs. 3 and 4. The labels "Side" and "Top" 
in the figures denote that the view from the direction normal and parallel to the rotating axis, 
respectively. Figure 3 shows the results for the relatively slow rotation with co = 0.025. The 
shape of the droplet viewed from the side becomes ellipse, while the shape is kept circle viewed 
from the top. This implies that the droplet deforms to an oblate spheroid. Figure 4 shows the 



Figure 3. Time evolutions of a droplet of = 16000 for u = 0.025. The yellow line presents 
the axis of the rotation. The views from the normal direction to the rotation axis are shown in 
the left, and those from the parallel direction are shown in the right. The form of the droplet 
changes from the sphere to the uniaxial. 




4. Same as Fig. 3 for u = 0.05. The form of the droplet becomes biaxial. 
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Figure 5. The angular velocity dependence of the deformation of rotating droplets. The radius 
ratio Rc/Rq is shown as a function of the reduced angular velocity oj/ujq. The solid line denotes 
the theoretical prediction in Eq. (25). 



results for the larger angular velocity u = 0.05. The both shapes viewed from side and top 
become ellipses, i.e., the droplet becomes biaxial shape. We find slight increase of temperature 
during rotation, but the increase is about ~ 1%. Therefore, the influence of the heating by 
rotation can be negligible. 

For the regime where the final shape of the droplet is uniaxial, we compare the results from 
MD and the theoretical predictions. We plot the ratio Rc/Rq as a function of the reduced 
frequency uj/ojq in Fig. 5. While the deviations from the theoretical prediction are observed, 
the difference between the numerical result and the theory becomes smaller as the size of the 
droplet increases. Note that, the experimental results such as in Ref. [8] shows good agreements 
with the theretical prediction, which implies that the size of the treated droplets is large enough 
so that the influence from the thickness of the surface is negligibly-small. 

4-2. Deformation Dynamics 

Next, we observe the dynamics of the deformation. While we have used the position of the Gibbs 
surface for the rotating droplet in steady state, it is difficult to determine it when the droplet 
undergoes deformation, since the density profile of such cases cannot be determined accurately. 
Therefore we measure the gyration radius instead of the position of the Gibbs surface in order 
to measure how a droplet deforms. The gyration radius R is defined as 

R'^^Url (28) 

where rj is the distance between particle i to the rotating axis. The time evolutions of the radius 
ratio R{t)/R[Q) for different values of oj are shown in Fig. 6. Deformation becomes larger for 
a larger value of uj. The time evolution of the radius is expected to be following exponential 
form as 

R{t) = i?eq - [i?eq " ^2(0)] exp(-t/r), (29) 
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Figure 6. Time evolutions of the deformation which is defined by the ratio of the current radius 
to the initial radius R{t)/R{0) for different values of the angular velocities. Temperature is fixed 
to be 0.6. 



with the final value of the radius Req, i.e., i?eq = linit-s.oo Rit). We define the reduced radius R* 

in order to investigate the characteristic timescale by removing the amplitude of the deformation. 
The time evolutions of the reduced radii are shown in Fig. 7. The data are collapsed into 
the single curve without scaling of the horizontal axis. This implies that the deformations 
for different rotation velocity have the same timescale. In other words, the timescale of the 
deformation is almost independent of the angular frequency. 

We also investigate infiuence of temperature. Time evolutions of the radii of the droplets for 
different temperatures are shown in Fig. 8. We estimate relaxation time r defined in Eq. (29), 
and determine the value of r = 88.0(3), 68.3(6), and 61.3(7) for temperatures T = 0.55, 0.60, 
and 0.70. One can see that the timescale becomes shorter as temperature increases. This means 
that the droplet deforms quickly for high temperature. 



5. Summary and Discussion 

To summarize, nonequilibrium deformation processes of rotating droplets are investigated by 
molecular dynamics simulations. First, we study the size-dependence of the deformations. While 
the deformations of the rotating droplets in steady state are smaller than those of the theoretical 
predictions, the differences decrease as the size of the droplets increases. The differences are the 
same order as the ratio of the thickness to the radius. Therefore, the differences come from the 
fact that the surface of a droplet has finite thickness, and the thickness cannot be ignored for 
small droplets. This difference was not observed in the past studies such as in Ref. [8], which 
implies that the treated droplet in the expriments are large enough. Note that, the finite-size 
effect does not comes from the number of particles consisting the droplet, but the ratio of length 
of the interface to the radius of the droplet. As the temperature increases, the length of the 
interface increases and diverges at the critical point. Therefore, the finite-size effect will be 
observed in the macroscopic experiments in the region near the critical point. For faster speed 
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Figure 7. Time evolutions of reduced radius R* for different values of the angular velocity u. 
They are well scaled. This implies that they share the identical timescale for their deformations, 
while their amplitudes of the deformations are different. 
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Figure 8. Time evolutions of deformations for different temperatures. The ratio of the radius 
R{t)/RQ is shown, where R{t) is the equatorial radius and R{0) is the initial radius, respectively, 
(inset) The temperature dependence of the relaxation time r. It shows that a droplet in higher 
temperature deforms quicker. 



of rotation, the uniaxial form becomes unstable and biaxial form is observed as reported in the 
past study using the finite-element method [6] . 

Nonequilibrium behaviors of the deformation are also observed. We find that the rotation 
frequency only affects the form of the droplet and does not change the timescale of deformations, 
when the rotation is small enough so that the final shapes are uniaxial. While the characteristic 
timescale is independent of the rotation speed, it decreases as temperature increases. In order 
to understand these nonequilibrium behaviors, we have constructed a following simple model. 
Consider a spring-mass system in fluid. A particle with mass M is connected to the spring 
whose spring constant is k, and the other end of the spring is fixed. When the particle moves 
with the angular velocity u around the fixed point, then an overdumped equation of motion is 
written as, 

T^ = M{xn + x)oj^ -kx, (31) 

where Xn is the natural length of the spring and x is displacement from it. F denotes a 
phenomenological dissipation coefficient which is proportional to viscosity of the fluid. The 
solution of the equation is 

x{t) = Xeq - [Xeq - x{0)] exp(-t/T) (32) 

with the length at the equilibrium 



and the relaxation time 



The solution Eq. (32) corresponds to the behavior Eq. (29). In this model, the spring constant 
k corresponds to the surface tension 7, which works as a restoring force (see Appendix). When 
the rotational frequency is much smaller than the restoring force, i.e., \uj'^\ <IC k/M, then the 
characteristic timescale is reduced to be, 

r^l (35) 

Then the timescale becomes independent of the rotational frequency, which is observed in the 
numerical simulations. The temperature dependence of the timescale r is not trivial, since 
both restoring force k and the dissipation coefficient T depend on temperature. Both surface 
tension and viscosity decreases as temperature increases, and the observed deformation timescale 
becomes shorter as temperature increases. This implies that the decrease of the viscosity is faster 
than that of the surface tension, while situations may change for other types of liquid. The 
dissipation coefficient F depends not only on temperature, but on size of a droplet. This size- 
dependence of dissipation is not trivial, since the dissipation caused by deformation is difficult to 
be solved exactly. These issues should be studied in the future. With larger angular velocity, the 
deformation of the droplet becomes more complicated. While rotation is almost rigid-body type 
for small angular velocity, the internal flow of droplets can change the behavior as discussed in 
Ref. [7]. Non rigid-body rotation will involve dissipation, which causes change of temperature, 
etc. Dynamics of such behaviors are also worth to be studied by molecular dynamics simulations. 
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Appendix A. Restoring Force and Surface Tension 

In this appendix, we relate the phenomenological restoring force k in Eq. (32) to the surface 
tension 7 and the rotating frequency a;. Consider a rotating droplet consisting particles. 
When the droplet rotate slow enough, the form of the droplet is uniaxial and well approximated 
by an oblate spheroid. The polar radius Rp and the equatorial radius R satisfy the following 
equation, 

47r 

N = -pR'Rp. (A.l) 

The surface area S of this droplet is 

/ „ i?^tanh~^e\ 
S = 2Tr{R^ + , (A.2) 



with the eccentricity e = ^1 — R^/R"^. When the rotation velocity is small enough, the 

eccentricity is very close to unity. Then the approximation tanh~^ e ~ e for |e| 1 leads 
to 

S = 27r(i?2 + Rl). (A.3) 
The moment of inertia / of the droplet is 

I = '^mNR^, (A.4) 

where m is the mass of the particles. Consequently, the effective potential of the droplet can be 
written as 

U = S-f-^Iuj^ (A.5) 

with the rotating frequency u and the surface tension 7. The contribution from the inner 
pressure is ignored here, since the observed pressure is almost constant during its rotation in 
the numerical results. The restoring force for the form of droplets F{r) is written as 

= - 47r7 ] R+ ' 
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The mechanical equilibrium radius is given by F(i?eq) = Oj which is 
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The linear approximation around the mechanical equilibrium point gives the effective spring 
constant of the droplet K as 



dR 



R—Rei 



247r7 . (A.8) 



This K corresponds to the restoring force k — Mu^ in Eq. (34). Then identifying mN with M, 
we have the characteristic timescale of this system as 



^ 247r7- 12MwV5' ^^'^^ 
When the rotation of the droplet is slow enough, the relaxation time becomes 

r = (A.IO) 

247r7 ^ ' 

which is essentially equivalent to Eq. (35). 
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